1.a) Se carga la informacion de los cajeros 108 al 116,
Caj108<-read.csv("cajero108.csv",header=F,dec=".",sep=";")
Caj109<-read.csv("cajero109.csv",header=F,dec=".",sep=";")
Caj110<-read.csv("cajero110.csv",header=F,dec=".",sep=";")
Caj111<-read.csv("cajero111.csv",header=F,dec=".",sep=";")
Caj112<-read.csv("cajero112.csv",header=F,dec=".",sep=";")
Caj113<-read.csv("cajero113.csv",header=F,dec=".",sep=";")
Caj114<-read.csv("cajero114.csv",header=F,dec=".",sep=";")
Caj115<-read.csv("cajero115.csv",header=F,dec=".",sep=";")
Caj116<-read.csv("cajero116.csv",header=F,dec=".",sep=";")
dim(Caj108)
## [1] 1 122
dim(Caj109)
## [1] 1 243
dim(Caj110)
## [1] 1 396
dim(Caj111)
## [1] 1 251
dim(Caj112)
## [1] 1 185
dim(Caj113)
## [1] 215 1
dim(Caj114)
## [1] 215 1
dim(Caj115)
## [1] 215 1
dim(Caj116)
## [1] 215 1
Convertimos ahora cada una en una serie de tiempo, debemos tener cuidado puesto que las cinco primeras se encuentran como vectores filas,
Caj108<-t(Caj108)
Caj109<-t(Caj109)
Caj110<-t(Caj110)
Caj111<-t(Caj111)
Caj112<-t(Caj112)
Caj108<-ts(Caj108[,1],start=c(2008,1),freq=365)
Caj109<-ts(Caj109[,1],start=c(2008,1),freq=365)
Caj110<-ts(Caj110[,1],start=c(2008,1),freq=365)
Caj111<-ts(Caj111[,1],start=c(2008,1),freq=365)
Caj112<-ts(Caj112[,1],start=c(2008,1),freq=365)
Caj113<-ts(Caj113[,1],start=c(2008,1),freq=365)
Caj114<-ts(Caj114[,1],start=c(2008,1),freq=365)
Caj115<-ts(Caj115[,1],start=c(2008,1),freq=365)
Caj116<-ts(Caj116[,1],start=c(2008,1),freq=365)
str(Caj108)
## Time-Series [1:122] from 2008 to 2008: 2440000 1248000 2499000 2079000 1335000 ...
## - attr(*, "names")= chr [1:122] "V1" "V2" "V3" "V4" ...
str(Caj109)
## Time-Series [1:243] from 2008 to 2009: 1989000 7780000 6859000 4802000 4733000 ...
## - attr(*, "names")= chr [1:243] "V1" "V2" "V3" "V4" ...
str(Caj110)
## Time-Series [1:396] from 2008 to 2009: 19222000 14495000 13983000 12565000 7124000 ...
## - attr(*, "names")= chr [1:396] "V1" "V2" "V3" "V4" ...
str(Caj111)
## Time-Series [1:251] from 2008 to 2009: 3562000 659000 2540000 1855000 1858000 ...
## - attr(*, "names")= chr [1:251] "V1" "V2" "V3" "V4" ...
str(Caj112)
## Time-Series [1:185] from 2008 to 2009: 4929000 1836000 1770000 1928000 4390000 ...
## - attr(*, "names")= chr [1:185] "V1" "V2" "V3" "V4" ...
str(Caj113)
## Time-Series [1:215] from 2008 to 2009: 24990 17780 18270 17260 20935 21810 9755 14620 10755 14000 ...
str(Caj114)
## Time-Series [1:215] from 2008 to 2009: 1470 265 855 435 735 255 600 170 715 430 ...
str(Caj115)
## Time-Series [1:215] from 2008 to 2009: 5450 3630 4710 3305 3590 4610 4630 2930 2470 3595 ...
str(Caj116)
## Time-Series [1:215] from 2008 to 2009: 1405 1810 885 890 2020 1615 1395 1850 1145 990 ...
Cargamos el paquete itsmr,
library(itsmr)
str(dowj)
## num [1:78] 111 111 110 111 111 ...
str(strikes)
## num [1:30] 4737 5117 5091 3468 4320 ...
str(Sunspots)
## num [1:100] 101 82 66 35 31 7 20 92 154 125 ...
str(wine)
## num [1:142] 464 675 703 887 1139 ...
Convertimos cada uno de los vectores de datos en series de tiempo. En este caso todos resultan ser vectores columnas y no exite informacion respecto a si las fechas se encuentran en listadas en orden creciente o decreciente por lo que se supondra que el orden de las fechas es exactamente el requerido
tdowj<-ts(dowj, start=c(1972, 241), freq=365)
plot(tdowj)
tstrikes<- ts(strikes, start=1951, end=1980, freq=1)
plot(tstrikes)
tSunspots<- ts(Sunspots, start=1870, end=1979, freq=1)
plot(tSunspots)
twine<- ts(wine, start=c(1980,1), end=c(1991,10), freq=12)
plot(twine)
str(tdowj)
## Time-Series [1:78] from 1973 to 1973: 111 111 110 111 111 ...
str(tstrikes)
## Time-Series [1:30] from 1951 to 1980: 4737 5117 5091 3468 4320 ...
str(tSunspots)
## Time-Series [1:110] from 1870 to 1979: 101 82 66 35 31 7 20 92 154 125 ...
str(twine)
## Time-Series [1:142] from 1980 to 1992: 464 675 703 887 1139 ...
Notemos que en el primer caso se utilizaron los par?metros \(1972\) y \(241\) para la fecha de iniciio de la serie, esto debido a que el 28 de agosto corresponde al dia 241 del a~no.
Cargamos la tabla de datos DJTable.csv,
DJTable<-read.csv("DJTable.csv",header=T,dec=".",sep=";")
tail(DJTable)
## X AA AXP BA BAC CAT CSCO CVX DD DIS GE
## 247 11/01/2010 17.45 41.47 60.87 16.93 64.13 24.59 80.88 34.26 31.36 16.76
## 248 08/01/2010 17.02 41.95 61.60 16.78 60.34 24.66 79.47 33.94 31.88 16.60
## 249 07/01/2010 16.61 41.98 62.20 16.93 59.67 24.53 79.33 34.39 31.83 16.25
## 250 06/01/2010 16.97 41.49 59.78 16.39 59.43 24.42 79.63 34.04 31.82 15.45
## 251 05/01/2010 16.13 40.83 58.02 16.20 59.25 24.58 79.62 33.93 31.99 15.53
## 252 04/01/2010 16.65 40.92 56.18 15.69 58.55 24.69 79.06 34.26 32.07 15.45
## HD HPQ IBM INTC JNJ JPM KO MCD MMM MRK MSFT
## 247 28.16 52.43 129.48 20.95 64.22 44.53 56.27 62.32 83.98 37.85 30.27
## 248 28.98 52.59 130.85 20.83 64.21 44.68 55.15 61.84 84.32 37.70 30.66
## 249 29.12 52.20 129.55 20.60 63.99 44.79 56.19 61.90 83.73 37.72 30.45
## 250 28.78 52.18 130.00 20.80 64.45 43.92 56.33 61.45 83.67 37.66 30.77
## 251 28.88 52.67 130.85 20.87 63.93 43.68 56.35 62.30 82.50 37.16 30.96
## 252 28.67 52.45 132.45 20.88 64.68 42.85 57.04 62.78 83.02 37.01 30.95
## PFE PG T TRV UNH UTX VZ WMT XOM
## 247 18.83 60.20 26.97 48.54 32.92 72.16 31.88 54.21 70.30
## 248 18.68 60.44 27.10 48.56 32.70 70.63 31.75 53.33 69.52
## 249 18.53 60.52 27.30 48.63 33.01 70.49 31.73 53.60 69.80
## 250 18.60 60.85 27.61 47.94 31.79 70.19 31.92 53.57 70.02
## 251 18.66 61.14 28.44 48.63 31.48 70.56 33.34 53.69 69.42
## 252 18.93 61.12 28.58 49.81 31.53 71.63 33.28 54.23 69.15
Vemos que las fechas estan en orden decreciente por lo que debemos invertirlas, para ello,
head(DJTable)
## X AA AXP BA BAC CAT CSCO CVX DD DIS GE
## 1 31/12/2010 15.39 42.92 65.26 13.34 93.66 20.23 91.25 49.88 37.51 18.29
## 2 30/12/2010 15.21 42.51 65.01 13.28 93.87 20.23 91.60 49.69 37.48 18.19
## 3 29/12/2010 15.13 42.86 65.05 13.31 93.78 20.25 91.37 50.02 37.60 18.27
## 4 28/12/2010 15.25 42.79 64.86 13.34 93.69 20.35 91.19 49.82 37.36 18.32
## 5 27/12/2010 15.23 43.05 64.75 13.27 94.07 20.16 90.12 49.63 37.48 18.19
## 6 23/12/2010 15.34 42.77 65.06 13.06 94.45 19.69 90.68 49.77 37.70 18.04
## HD HPQ IBM INTC JNJ JPM KO MCD MMM MRK MSFT PFE
## 1 35.06 42.10 146.76 21.03 61.85 42.42 65.77 76.76 86.30 36.04 27.91 17.51
## 2 34.86 42.26 146.67 21.02 61.94 42.23 65.50 76.76 86.54 36.01 27.85 17.49
## 3 34.89 42.32 146.52 20.94 62.13 42.36 65.45 76.99 86.76 36.21 27.97 17.60
## 4 35.09 42.25 145.71 20.88 62.05 42.61 65.36 76.43 86.74 36.20 28.01 17.59
## 5 35.24 41.82 145.34 20.84 61.93 42.67 65.07 76.43 87.01 36.23 28.07 17.49
## 6 35.09 41.74 145.89 20.84 62.25 42.08 65.58 76.96 86.47 36.29 28.30 17.61
## PG T TRV UNH UTX VZ WMT XOM
## 1 64.33 29.38 55.71 36.11 78.72 35.78 53.93 73.12
## 2 64.28 29.33 55.54 35.94 78.85 35.56 54.07 73.36
## 3 64.40 29.31 55.59 35.91 79.10 35.58 54.08 73.37
## 4 64.76 29.23 55.66 35.72 79.29 35.62 53.74 73.42
## 5 64.67 29.25 55.79 35.54 79.27 35.50 53.57 73.01
## 6 65.24 29.20 55.48 35.77 79.50 35.44 53.60 73.20
CSCO<-rev(DJTable$CSCO)
CVX<-rev(DJTable$CVX)
DD<-rev(DJTable$DD)
Por ultimo las convertimos en series de tiempo,
CSCO<-ts(CSCO, start=c(2010, 1), freq=365)
CVX<-ts(CVX, start=c(2010, 1), freq=365)
DD<-ts(DD, start=c(2010, 1), freq=365)
str(CSCO)
## Time-Series [1:252] from 2010 to 2011: 24.7 24.6 24.4 24.5 24.7 ...
str(CVX)
## Time-Series [1:252] from 2010 to 2011: 79.1 79.6 79.6 79.3 79.5 ...
str(DD)
## Time-Series [1:252] from 2010 to 2011: 34.3 33.9 34 34.4 33.9 ...
Cargamos la informacion del cajero 103,
Caj103<-read.csv("cajero103.csv",header=F,dec=".",sep=";")
dim(Caj103)
## [1] 1065 5
Debemos transformar a Caj103 en un vector columna, para ello
Caj103<-as.matrix(Caj103) ### Pues Caj103 es un DataFrame
Caj103<-as.vector(Caj103)
length(Caj103)
## [1] 5325
Convertimos ahora el vector en una serie de tiempo
Caj103<-ts(Caj103,start=c(1998,1),freq=365)
plot(Caj103)
Para verificar la normalidad graficamos el histograma para los residuos,
hist(diff(Caj103),prob=T,ylim=c(0,2.5e-07),col="red")
lines(density(diff(Caj103)),lwd=2)
mu<-mean(diff(Caj103))
sigma<-sd(diff(Caj103))
x<-seq(-1e+07,1e+07,length=1000)
y<-dnorm(x,mu,sigma)
lines(x,y,lwd=2,col="blue")
Vemos de la grafica que el supuesto de normalidad se verifica aproximadamente bien.
Realizamos el suavizamiento de la serie para \(a=4\), \(a=6\) y \(a=10\),
plot(Caj103, type="l")
st.1<- filter(Caj103,filter=rep(1/9,9))
st.2<-filter(Caj103,filter=rep(1/13,13))
st.3<-filter(Caj103,filter=rep(1/21,21))
lines(st.1,col="red")
lines(st.2,col="purple")
lines(st.3,col="blue")
Esta serie presenta mucha variabilidad, vemos como ella se reduce al hacer el suavizamiento o filtrado.
Realizamos la dscomposicion de la serie utilizando la funcion stl. Para ello utilizamos la serie de tiempo Caj103 que ya ha sido generada,
DescCaj103<-stl(Caj103,s.window="periodic")
head(DescCaj103$time.series)
## seasonal trend remainder
## [1,] 1831668.00 3350290 -4553958.1
## [2,] 809017.04 3351466 -2998482.8
## [3,] 359299.40 3352641 179059.2
## [4,] 72981.77 3353817 -421798.8
## [5,] -65002.53 3354993 -573990.2
## [6,] -571453.49 3356168 -163714.9
plot(DescCaj103)
Cargamos el paquete itsmr y cargamos la tabla de datos deaths,
suppressMessages(library(itsmr))
str(deaths)
## num [1:72] 9007 8106 8928 9137 10017 ...
Transformamos ahora la tabla en una serie de tiempo. Es importante recordar que los datos corresponden a muertes accidentales en Estados Unidos entre los a?os \(1973\) y \(1978\). Dado que la tabla tiene \(72\) entradas la frecuencia de losd datos es mensual, por lo tanto la frecuencia es igual a \(12\).
Deaths<-ts(deaths, start=c(1973,1), freq=12)
Ahora bien, calculamos primero el peridiograma
res<-spec.pgram(Deaths, log = "no")
res.ord<-order(res$spec,res$freq,decreasing=TRUE)
res.ord[1:4]
## [1] 6 12 1 30
Vemos que los m?ximos se encuentran en las posiciones \(6, 12 y 30\).
max1<-res$freq[6]
period1<-12/max1
max2<-res$freq[12]
period2<-12/max2
max3<-res$freq[30]
period3<-12/30
c(period1,period2,period3)
## [1] 12.0 6.0 0.4
res<-spec.pgram(Deaths, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")
Vemos que el primer periodo es \(12\), esto justifica el utilizar el modelo SARIMA con perdiodo \(12\).
Generamos el modelo y hacmeos la predicci?n. No quedaba claro del enunciado si una predicci?n correspond?a a predecir s?lo un nuevo dato o a realizar una predicci?n a voluntad. En este caso se eligi? predecir un a?o,
fit<-arima(Deaths,order=c(1,2,1),seasonal=list(order=c(2,1,2),period=12))
LH.pred<-predict(fit,n.ahead=12)
layout(matrix(c(1,1,2,3),2,2,byrow=TRUE))
plot(Deaths,xlim=c(1973,1980),ylim=c(7000,12000),type="o")
lines(LH.pred$pred,col="green",type="o")
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o")
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Deaths,main="Autocorrelaci?n?Simple",col="black",ylim=c(-1,1))
pacf(Deaths,main="Autocorrelaci?n?Parcial",col="black",ylim=c(-1,1))
Cargamos la tabla de datos correspondiente al cajero \(101\). Recordemos que debemos trasponer este vector de datos para luego transformarlo en una serie de tiempo. La serie comienza el \(1\) de enero del a?o \(2008\), y tiene frecuencia diaria.
Caj101<-read.csv("Cajero101.csv",header=F,dec=".",sep=";")
Caj101<-t(Caj101)#?Matriz?transpuesta
Caj101<-ts(Caj101[,1],start=c(2008,1),freq=365)
plot(Caj101,type="o",col="blue")
Ahora bien, calculamos primero el peridiograma
res<-spec.pgram(Caj101, log = "no")
res.ord<-order(res$spec,res$freq,decreasing = TRUE)
res.ord[1:4]
## [1] 101 439 220 438
Vemos que los m?ximos se encuentran en las posiciones \(101, 439\) y \(220\).
max1<-res$freq[101]
period1<-365/max1
max2<-res$freq[439]
period2<-365/max2
max3<-res$freq[220]
period3<-365/max3
c(period1,period2,period3)
## [1] 15.207921 3.498861 6.981818
res<-spec.pgram(Caj101, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")
Vemos que el primer periodo es \(15\), este es el periodo m?s adecuado. Realizamos la predicci?n ajustando un modelo \(SARIMA(1,1,2)(2,1,2)\) con per?odo \(15\).
fit<-arima(Caj101,order=c(1,1,2),seasonal=list(order=c(2,1,2),period=15))
LH.pred<-predict(fit,n.ahead=365)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(Caj101,xlim=c(2008,2013),ylim=c(min(Caj101)-3500000,max(Caj101)+1000000),type="o")
lines(LH.pred$pred,col="green",type="o")
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o")
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Caj101, main="Autocorrelaci?n Simple",col="black",ylim=c(-1,1))
pacf(Caj101,main="Autocorrelaci?n Parcial",col="black",ylim=c(-1,1))
Ahora bien, extraemos los ?ltimos dos meses de la serie,
Caj101<-Caj101[1461:1521]
Caj101<-ts(Caj101,start=c(2012,2),freq=365)
plot(Caj101,type="o",col="blue")
Ahora bien, calculamos el peridiograma
res<-spec.pgram(Caj101, log = "no")
res.ord<-order(res$spec,res$freq,decreasing = TRUE)
res.ord[1:4]
## [1] 2 18 19 1
Vemos que los m?ximos se encuentran en las posiciones \(2, 18\) y \(19\).
max1<-res$freq[2]
period1<-365/max1
max2<-res$freq[18]
period2<-365/max2
max3<-res$freq[19]
period3<-365/max3
c(period1,period2,period3)
## [1] 32.000000 3.555556 3.368421
res<-spec.pgram(Caj101, log = "no")
abline(v=max1, lty="dotted",col="red")
abline(v=max2, lty="dotted",col="blue")
abline(v=max3, lty="dotted",col="magenta")
Vemos que el primer periodo es \(32\), pero este periodo es demasiado grande dado que los datos corresponden s?lo a dos meses. Por esto se elegir? el siguiente periodo que es iguala \(3\).
fit<-arima(Caj101,order=c(7,1,7),seasonal=list(order=c(2,1,2),period=3))
## Warning in arima(Caj101, order = c(7, 1, 7), seasonal = list(order = c(2, :
## possible convergence problem: optim gave code = 1
LH.pred<-predict(fit,n.ahead=30)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(Caj101,xlim=c(2012,2012.25),ylim=c(min(Caj101)-5000000,max(Caj101)+10000000),type="o")
lines(LH.pred$pred,col="green",type="o")
lines(LH.pred$pred+2*LH.pred$se,col="green",lty=3,type="o")
lines(LH.pred$pred-2*LH.pred$se,col="green",lty=3,type="o")
acf(Caj101, main="Autocorrelaci?n Simple",col="black",ylim=c(-1,1))
pacf(Caj101,main="Autocorrelaci?n Parcial",col="black",ylim=c(-1,1))
Vemos que en este caso la prediccion parece tener un comportamiento un mas dinamico y no solo a describir la media.
Comenzaremos por generar las funciones que nos permitir?n calcular los diferentes tipos de error y calibrar el modelo de Holt-Winters.
La siguiente funci?n nos permitir? calcular diferentes tipos de error, Error Relativo, Error Cuadr?tico Medio, Total de Predicciones Sobre el Valor Real, Error Relativo para las Predicciones por Sobre el Valor Real.
La siguiente funci?n nos permitir? calibrar el modelo de Holt-Winters. Se modific? la funci?n presentada en el curso de modo de elegir el incremento como par?metro de la funci?n.
calibrar<-function(serie.aprendizaje,serie.testing, pot.incremento) {
incremento<-10^{pot.incremento}
error.c<-Inf
alpha.i<-incremento # alpha no puede ser cero
while(alpha.i<=1) {
beta.i<-0
while(beta.i<=1) {
gamma.i<-0
while(gamma.i<=1) {
mod.i<-HoltWinters(serie.aprendizaje,alpha=alpha.i,beta=beta.i,gamma=gamma.i)
res.i<-predict(mod.i,n.ahead=length(serie.testing))
error.i<-sqrt(ECM(res.i,serie.testing))
if(error.i<error.c) {
error.c<-error.i
mod.c<-mod.i
}
gamma.i<-gamma.i+incremento
}
beta.i<-beta.i+incremento
}
alpha.i<-alpha.i+incremento
}
return(mod.c)
}
Cargamos el paquete itsmr y cargamos la tabla de datos deaths,
suppressMessages(library(itsmr))
str(deaths)
## num [1:72] 9007 8106 8928 9137 10017 ...
Transformamos ahora la tabla en una serie de tiempo. Es importante recordar que los datos corresponden a muertes accidentales en Estados Unidos entre los a?os \(1973\) y \(1978\). Dado que la tabla tiene \(72\) entradas la frecuencia de losd datos es mensual, por lo tanto la frecuencia es igual a \(12\).
Deaths<-ts(deaths, start=c(1973,1), freq=12)
Construimos las tablas de aprendizaje y prueba, donde esta ?ltima contiene los ?ltimos diez datos:
Deaths.train<-Deaths[1:62]
Deaths.train<-ts(Deaths.train, start=c(1973,1), freq=12)
Deaths.test<-Deaths[63:72]
Deaths.test<-ts(Deaths.test, start=c(1978,3), freq=12)
Calibramos el moldelo de Holt-Winters para la serie de datos de aprendizaje:
modelo.HW<-calibrar(Deaths.train,Deaths.test, -1)
modelo.HW
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = serie.aprendizaje, alpha = alpha.i, beta = beta.i, gamma = gamma.i)
##
## Smoothing parameters:
## alpha: 0.1
## beta : 0.3
## gamma: 0.2
##
## Coefficients:
## [,1]
## a 8607.57551
## b 26.30563
## s1 -710.30464
## s2 -407.20322
## s3 216.10027
## s4 792.95226
## s5 1750.51793
## s6 1076.09047
## s7 55.89459
## s8 441.49007
## s9 -111.42316
## s10 47.73263
## s11 -876.26319
## s12 -1617.78781
Ahora bien, recordamos de la tarea anterior que el periodo asociado a estos datos es 12. Calibrando el modelo ARIMA mediante el uso de la funci?n \(auto.arima\) obtenemos:
suppressMessages(library(forecast))
auto.arima(Deaths.train)
## Series: Deaths.train
## ARIMA(0,1,1)(0,1,1)[12]
##
## Coefficients:
## ma1 sma1
## -0.4295 -0.4787
## s.e. 0.1368 0.2061
##
## sigma^2 estimated as 116720: log likelihood=-356.01
## AIC=718.03 AICc=718.56 BIC=723.7
Donde vemos que los par?metros ?ptimos a utilizar para ajustar el modelo ARIMA son \((0,1,1)(0,1,1)\), construyendo el modelo,
modelo.A<-arima(Deaths.train,order=c(0,1,1),seasonal=list(order=c(0,1,1),period=12))
Generamos \(10\) meses de predicciones utilizando ambos modelos. Primero, la predicci?n correspondiente al m?todo de Holt-Winters,
predict.HW<-predict(modelo.HW,n.ahead=10)
plot(predict.HW,xlim=c(1973,1979), ylim=c(6800, 11500), main = "Predicciones mediante modelo de Holt-Winters", col=2)
lines(Deaths.train,col=1)
y la predicci?n correspondiente al modelo SARIMA deperiodo \(12\) y par?metros \(()()\)
predict.A<-predict(modelo.A,n.ahead=10)
plot(predict.A$pred,xlim=c(1973,1979), ylim=c(6800, 11500), main = "Predicciones mediante modelo SARIMA", col = 2)
lines(Deaths.train,col=1)
Calculamos ahora el Error Cuadr?tico Medio asociado a cada predicci?n:
ECM.HW<-sqrt(ECM(predict.HW,Deaths.test))
ECM.A<-sqrt(ECM(predict.A$pred,Deaths.test))
Errores.ECM<-c(ECM.HW, ECM.A)
names(Errores.ECM)<-c("Holt-Winters", "SARIMA")
Errores.ECM
## Holt-Winters SARIMA
## 173.6803 494.5145
Vemos claramente que el modelo de Holt-Winters es mejor que el SARIMA, esto puesto que fue calibrado precisamente para minimizar el error cuadr?tico medio. Si utiliz?ramos otra medida de error podr?a darse una situaci?n diferente. De todas formas, la misma funci?n podria utilizarse para minimizar otros tipos de error.
Ahora bien, el Error Relativo para ambos modelos est? dado por:
ER.HW<-ER(predict.HW,Deaths.test)
ER.A<-ER(predict.A$pred,Deaths.test)
Errores.R<-c(ER.HW, ER.A)
names(Errores.R)<-c("Holt-Winters", "SARIMA")
Errores.R
## Holt-Winters SARIMA
## 0.01579028 0.04676109
En este caso observamos tambi?n que es el modelo de Holt-Winters el que tiene asociado un menor error. El mejor modelo basado en este criterio es entonces el modelo de Holt-Winters.
Cargamos la tabla de datos correspondiente al cajero \(101\). Recordemos que debemos trasponer este vector de datos para luego transformarlo en una serie de tiempo. La serie comienza el \(1\) de enero del a?o \(2008\), y tiene frecuencia diaria.
Caj101<-read.csv("Cajero101.csv",header=F,dec=".",sep=";")
Caj101<-t(Caj101)#?Matriz?transpuesta
Caj101<-ts(Caj101[,1],start=c(2008,1),freq=365)
plot(Caj101,type="o",col="blue")
str(Caj101)
## Time-Series [1:1521] from 2008 to 2012: 5409000 4703000 4243000 6395000 4645000 ...
## - attr(*, "names")= chr [1:1521] "V1" "V2" "V3" "V4" ...
Repetimos los pasos seguidos en el ejercicio anterior, esto es, construimos las tablas de aprendizaje y prueba, donde esta ?ltima contiene los ?ltimos 31 d?as:
Caj101.train<-Caj101[1:1490]
Caj101.train<-ts(Caj101.train, start=c(2008,1), freq=365)
Caj101.test<-Caj101[1491:1521]
Caj101.test<-ts(Caj101.test, start=c(2012,30), freq=365)
Calibramos el moldelo de Holt-Winters para la serie de datos de aprendizaje:
modelo.HW<-calibrar(Caj101.train,Caj101.test, -1)
modelo.HW
## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = serie.aprendizaje, alpha = alpha.i, beta = beta.i, gamma = gamma.i)
##
## Smoothing parameters:
## alpha: 0.1
## beta : 0
## gamma: 0.5
##
## Coefficients:
## [,1]
## a 2271522.4035
## b 854.2779
## s1 -57574.8707
## s2 1530778.3887
## s3 314372.0384
## s4 350133.6202
## s5 434520.2718
## s6 -866898.5304
## s7 -1155317.7307
## s8 -484213.3089
## s9 280161.0677
## s10 -438313.8717
## s11 59972.5202
## s12 356671.6813
## s13 -1046435.2936
## s14 -744261.8615
## s15 1897087.0799
## s16 2310506.3525
## s17 198026.1161
## s18 -13928.9809
## s19 335633.4082
## s20 -732956.9313
## s21 -725233.2713
## s22 -410112.8102
## s23 -593490.9936
## s24 -1077907.6844
## s25 1175.2000
## s26 224114.4342
## s27 -522736.1947
## s28 -1194135.3761
## s29 -67387.9964
## s30 1186733.6711
## s31 1836692.9823
## s32 1030195.6703
## s33 1371681.5556
## s34 -252439.3857
## s35 -87414.9217
## s36 611178.6592
## s37 -195914.3419
## s38 -44589.1168
## s39 -63857.2025
## s40 -364342.5538
## s41 -1270778.2540
## s42 -1187134.2767
## s43 -685054.7330
## s44 698706.8585
## s45 383887.6042
## s46 1937734.5588
## s47 492766.5719
## s48 215962.4992
## s49 104531.3040
## s50 -77266.6927
## s51 -303137.5129
## s52 -579220.5384
## s53 -1239284.8362
## s54 -992544.0477
## s55 -1417499.8884
## s56 -1436511.5520
## s57 -823162.3830
## s58 20380.7885
## s59 3145156.1260
## s60 880246.2415
## s61 -473229.9273
## s62 919115.1425
## s63 -522243.8008
## s64 -672140.2023
## s65 -755060.4137
## s66 -471695.2599
## s67 -553080.1691
## s68 28628.6654
## s69 -1366367.6798
## s70 -885176.8804
## s71 -501991.8133
## s72 1326950.6865
## s73 89627.3050
## s74 1248168.1676
## s75 400711.5532
## s76 -843084.6366
## s77 640771.2888
## s78 137818.7581
## s79 64853.2447
## s80 -427491.5093
## s81 118539.4466
## s82 -431168.4284
## s83 -1537357.0308
## s84 -1275983.2206
## s85 -657962.6687
## s86 506166.1566
## s87 -1172892.1328
## s88 -1463692.8849
## s89 359786.0696
## s90 -610899.8290
## s91 -396678.7524
## s92 779821.2484
## s93 1155616.3219
## s94 -125192.6910
## s95 532472.1697
## s96 318850.0210
## s97 -1204803.0431
## s98 -1284592.0140
## s99 -818018.3294
## s100 760703.0856
## s101 -1650923.1224
## s102 -334419.0467
## s103 -1317132.8392
## s104 -2205939.0952
## s105 -214622.3190
## s106 905171.3820
## s107 440205.9355
## s108 -268751.6520
## s109 -861325.6783
## s110 -256006.2746
## s111 -1137266.4437
## s112 -1164405.8828
## s113 -1668804.4763
## s114 -1960150.9449
## s115 -1590926.7648
## s116 -1024037.9947
## s117 -32585.5337
## s118 -1596228.8421
## s119 -1218303.3608
## s120 1219534.8313
## s121 1050806.6307
## s122 1424433.0484
## s123 -988827.9473
## s124 -688577.0783
## s125 -895305.7415
## s126 -916996.6005
## s127 -281295.1088
## s128 -552943.6189
## s129 -576499.3522
## s130 -772528.9856
## s131 726592.9211
## s132 -1495977.2418
## s133 -562624.2538
## s134 552914.7338
## s135 1538591.6907
## s136 -159049.8362
## s137 -283167.8556
## s138 1339665.5038
## s139 -2061466.3671
## s140 -928108.0471
## s141 -523478.1660
## s142 -896653.6772
## s143 -980748.3795
## s144 -513257.5095
## s145 -339707.8401
## s146 -1744099.6188
## s147 -1519503.0303
## s148 -820309.5803
## s149 -263554.0469
## s150 537116.8739
## s151 1136530.7701
## s152 -51521.1827
## s153 -1364749.6578
## s154 -226667.5328
## s155 50031.6830
## s156 -493198.7546
## s157 -978270.1957
## s158 -101482.7305
## s159 -838761.1843
## s160 -1838976.7325
## s161 -1559949.4432
## s162 -908167.5582
## s163 -1401897.5196
## s164 379090.3464
## s165 928976.4844
## s166 1208967.0285
## s167 -945528.8514
## s168 459766.5331
## s169 608085.0186
## s170 -250707.6197
## s171 -723705.7659
## s172 -786469.0734
## s173 -634019.8944
## s174 -1221759.8176
## s175 -1483245.2442
## s176 -786543.2990
## s177 -984328.6882
## s178 -161510.2295
## s179 2019913.0496
## s180 1192017.7670
## s181 -297780.9912
## s182 -199565.5525
## s183 431219.5071
## s184 848933.3095
## s185 391177.1450
## s186 284714.4674
## s187 -130446.8790
## s188 -1200497.9256
## s189 -1247719.8615
## s190 -639552.0267
## s191 -739456.2777
## s192 -648597.7593
## s193 856233.7560
## s194 1165767.3946
## s195 -1026072.8610
## s196 143917.9949
## s197 100030.3470
## s198 1050055.5615
## s199 292348.0507
## s200 210470.9653
## s201 -402832.5320
## s202 -1054224.8482
## s203 -1190383.4195
## s204 -599996.1032
## s205 -1156013.5035
## s206 -1161336.1664
## s207 159195.0016
## s208 -349217.0999
## s209 -1551269.3683
## s210 -1170059.4173
## s211 525929.8008
## s212 1458934.4663
## s213 1724563.9160
## s214 -296335.1809
## s215 496886.7791
## s216 -746656.1970
## s217 2107.5115
## s218 -603547.2638
## s219 -385039.7857
## s220 -15107.2559
## s221 28646.1263
## s222 -554288.5892
## s223 -1435195.2973
## s224 -974264.4597
## s225 721813.4350
## s226 1690289.4023
## s227 1213113.6627
## s228 892992.6557
## s229 931711.4010
## s230 -721573.3761
## s231 200583.5348
## s232 -689677.6826
## s233 -214446.4390
## s234 -409734.4339
## s235 -446508.7964
## s236 -833243.4596
## s237 -1325624.7615
## s238 -1112237.0203
## s239 -1094700.6780
## s240 125554.9604
## s241 -362367.6669
## s242 1024495.6129
## s243 894839.1808
## s244 -1058108.4873
## s245 1006801.0201
## s246 1340146.5082
## s247 299666.3034
## s248 31225.3069
## s249 -51598.7067
## s250 -392685.8589
## s251 -1015582.6467
## s252 -1309427.7494
## s253 -813623.7402
## s254 -799871.0730
## s255 -391218.0587
## s256 1981980.1301
## s257 96744.7524
## s258 -1014427.6722
## s259 1599283.3060
## s260 844586.5270
## s261 248962.3226
## s262 -46900.6243
## s263 -258602.0729
## s264 -421730.2924
## s265 -951342.6386
## s266 -1042563.7235
## s267 -705933.0517
## s268 -761829.8545
## s269 -376244.7857
## s270 957613.0287
## s271 316849.6330
## s272 -878773.1371
## s273 680437.0065
## s274 645861.3716
## s275 995660.5598
## s276 1123135.7019
## s277 182811.5573
## s278 306185.4402
## s279 -1237698.3390
## s280 -863474.0847
## s281 -1122166.5140
## s282 -84173.3579
## s283 129642.1416
## s284 -364744.4754
## s285 -879511.8333
## s286 -1451715.1767
## s287 -156703.6883
## s288 1497951.1736
## s289 952692.4533
## s290 520701.5986
## s291 442103.6895
## s292 -1016546.4413
## s293 -1382460.3225
## s294 -734269.7635
## s295 -1163797.3848
## s296 -886136.7049
## s297 -826968.1253
## s298 -214610.7451
## s299 -447057.1887
## s300 -857719.3093
## s301 -423112.3321
## s302 -746496.8626
## s303 679345.0459
## s304 2794777.9016
## s305 2390695.0981
## s306 -149021.9089
## s307 118255.8809
## s308 26179.2061
## s309 -235481.6333
## s310 -375776.1404
## s311 -369359.2738
## s312 -426444.4102
## s313 -822409.5428
## s314 -1215928.3229
## s315 -1135169.6542
## s316 -631545.0012
## s317 453019.9995
## s318 1171214.2607
## s319 2753812.4202
## s320 -147868.5510
## s321 402593.8243
## s322 409111.0680
## s323 -232616.2614
## s324 -347333.8329
## s325 218011.7302
## s326 136338.5645
## s327 -311262.4753
## s328 -957009.3846
## s329 -345387.8165
## s330 -205125.5307
## s331 624859.5294
## s332 1579114.5696
## s333 2580632.4006
## s334 1747073.7404
## s335 785464.2236
## s336 2661759.3267
## s337 2669837.2039
## s338 1443485.1093
## s339 1278156.4800
## s340 2413080.7979
## s341 799938.1078
## s342 1288736.6400
## s343 1873821.8572
## s344 1162773.3790
## s345 1875684.0414
## s346 2325187.2224
## s347 3053383.1894
## s348 1998182.9231
## s349 2217575.3548
## s350 1668473.1205
## s351 2638687.4137
## s352 2669125.6905
## s353 1629548.1397
## s354 2632874.1944
## s355 1223346.7772
## s356 1236551.7183
## s357 1423710.8372
## s358 1548086.6157
## s359 -932974.6608
## s360 983174.3889
## s361 685792.3608
## s362 313177.5169
## s363 -642116.0335
## s364 413025.6429
## s365 258159.7706
Ahora bien, recordamos de la tarea anterior que el periodo asociado a estos datos es \(15\). Calibrando el modelo ARIMA mediante el uso de la funci?n \(auto.arima\) obtenemos:
suppressMessages(library(forecast))
auto.arima(Caj101.train)
## Series: Caj101.train
## ARIMA(2,1,2)
##
## Coefficients:
## ar1 ar2 ma1 ma2
## -0.4580 0.1426 -0.1807 -0.7834
## s.e. 0.0899 0.0483 0.0832 0.0823
##
## sigma^2 estimated as 2.838e+12: log likelihood=-23460.14
## AIC=46930.27 AICc=46930.31 BIC=46956.8
Donde vemos que los par?metros ?ptimos a utilizar para ajustar el modelo ARIMA son \((2,1,2)(2,1,2)\), construyendo el modelo,
modelo.A<-arima(Caj101.train,order=c(2,1,2),seasonal=list(order=c(2,1,2),period=15))
Generamos \(31\) d?as de predicciones utilizando ambos modelos. Primero, la predicci?n correspondiente al m?todo de Holt-Winters,
predict.HW<-predict(modelo.HW,n.ahead=31)
plot(Caj101.train,xlim=c(2008,2012.1), main = "Predicciones mediante modelo de Holt-Winters")
lines(predict.HW,col=2)
y la predicci?n correspondiente al modelo SARIMA de periodo \(15\) y par?metros \((2,1,2)(2,1,2)\)
predict.A<-predict(modelo.A,n.ahead=10)
plot(Caj101.train,xlim=c(2008,2012.1), main = "Predicciones mediante modelo SARIMA", col = 1)
lines(predict.A$pred,col=2)
Calculamos ahora el Error Cuadr?tico Medio asociado a cada predicci?n:
ECM.HW<-sqrt(ECM(predict.HW,Caj101.test))
ECM.A<-sqrt(ECM(predict.A$pred,Caj101.test))
Errores.ECM<-c(ECM.HW, ECM.A)
names(Errores.ECM)<-c("Holt-Winters", "SARIMA")
Errores.ECM
## Holt-Winters SARIMA
## 1501120.4 739513.6
Vemos claramente que el modelo SARIMA tiene un error cuadr?tico medio m?s bajo que el de Holt-Winters.
Ahora bien, el Error Relativo para ambos modelos est? dado por:
ER.HW<-ER(predict.HW,Caj101.test)
ER.A<-ER(predict.A$pred,Caj101.test)
Errores.R<-c(ER.HW, ER.A)
names(Errores.R)<-c("Holt-Winters", "SARIMA")
Errores.R
## Holt-Winters SARIMA
## 0.4070711 0.1319049
En este caso observamos que, al igual que en el caso del error cuadr?tico medio, el error relativo asociado al modelo SARIMA es menor que el asociado al m?todo de Holt-Winters
Basados en ambos errores, el modelo SARIMA parece ser el m?s adecuado para modelar estos datos. Se elije entonces este modelo para realizar la predicci?n para siete d?as. Debemos primero generar el modelo para todo el conjunto de datos,
modelo.A<-arima(Caj101,order=c(2,1,2),seasonal=list(order=c(2,1,2),period=15))
Ahora realizamos la predicci?n para siete d?as,
Caj101.pred<-predict(modelo.A,n.ahead=7)
layout(matrix(c(1,1,2,3), 2, 2, byrow = TRUE))
plot(Caj101,xlim=c(2008,2012.1),ylim=c(0, max(Caj101)+1000000),type="o")
lines(Caj101.pred$pred,col=2)
Una idea, a priori, respecto a la calidad de esta predicci?n puede generarse gracias al c?lculo de los tipos de error asociados al modelo. Este c?lculo no utiliza los valores reales para los siete d?as predichos y por ende no es una medida necesariamente precisa. A posteriori, esto es, una vez conocido el valor de la serie para algunos de estos siete d?as, podr?a calcularse alg?n estimado del error para evaluar la predicci?n en base a esta nueva informaci?n y con ello decidir si esta predicci?n ha sido pertinente o no. Recordemos que lo ?ptimo seria actualizar las predicciones cada vez que contemos con nueva informaci?n puesto que mientras m?s datos sean predecidos, mayor es la varianza asociada a esta predicci?n y por ende mayor probabilidad de que el valor predicho se encuentre lejos del valor real.
Cargamos la tabla de datos del Dow Jones (DJTable.csv), y la tabla traspuesta (DJTableTranspose.csv),
DJTable<-read.csv("DJTable.csv",header=T,dec=".",sep=";")
DJTableTranspose<-read.csv("DJTableTranspose.csv",header=T,dec=".",sep=";")
Realizamos un ACP para la tabla DJTable,
suppressMessages(library(FactoMineR))
res<-PCA(DJTable[,-1], scale.unit=TRUE, ncp=6, graph = FALSE)
plot(res, axes=c(1, 2), choix="ind", col.ind="Red",new.plot=TRUE)
plot(res, axes=c(1, 2), choix="var", col.var="blue",new.plot=TRUE)
y para la tabla traspuesta,
res<-PCA(DJTableTranspose[,-1], scale.unit=TRUE, ncp=6, graph = FALSE)
plot(res, axes=c(1, 2), choix="ind", col.ind="Red",new.plot=TRUE)
plot(res, axes=c(1, 2), choix="var", col.var="blue",new.plot=TRUE)
En el primer caso detectamos al menos tres grupos diferentes. Mirando el c?rculo de correlaciones vemos que cada uno corresponde a combinaciones de indicadores que tomaron valores altos ese d?a. Dado que en este caso los datos corresponden a valores de \(30\) indocadores econ?micos, no es sencillo realizar la interpretaci?n puesto que aparecen muchos efectos posiblemente mezclados.
En el segundo caso, se ha realizado un ACP tomando como atributos a la fecha en la que los correspondiente indicadores econ?micos tomaron los correspondientes valores. En este caso, pueden identificarse dos marcados grupos, el compuesto por el indicador \(IBM\) y otro compuesto por los dem?s. Analizando el c?rculo de correlaciones, concluimos que el primer grupo est? conformado por un indicador que toma valores altos en todos los d?as.
Recordemos que los datos ya fueron cargados en memoria durante la resoluci?n del ejercico anterior. En este caso, realizamos clasificaci?n jer?rquica utilizando el paquete dtw,
suppressWarnings(suppressMessages(library(dtw)))
suppressWarnings(suppressMessages(library(rattle)))
clust.series <- hclust(dist(DJTable[,-1],method="DTW"), method="average")
par(mfrow=c(1,1))
plot(clust.series)
Graficamos ahora los centros de gravedad para los tres grupos,
centros<-centers.hclust(DJTable[,-1],clust.series,nclust=3,use.median=FALSE)
par(mfrow=c(3,1))
plot(centros[1,],type="o")
plot(centros[2,],type="o")
plot(centros[3,],type="o")
Las gr?ficas parecen ser bastante similares, dado que la escala de las gr?ficas es diferente, se mostrar?n en el mismo plot de modo de poder compararlas de forma m?s objetiva
par(mfrow=c(1,1))
plot(centros[1,],type="o")
lines(centros[2,],type="o", col=2)
lines(centros[3,],type="o", col=3)
Vemos que el primer cl?ster est? formado por aquellas fechas en las que los indicadores \(CAT\), \(CVX\), \(IBM\), \(KO\), \(MCD\) y \(MMM\) tomaron valores considerablemente m?s altos. El segundo cl?ster est? formado por aquellas fechas en las que los indicadores \(AA\), \(PFE\), \(UNH\), \(VZ\), \(WMT\) y \(XOM\) alcanzaron los valores m?s bajos, por ?ltimo, el tercer cl?ster est? formado por aquellas fechas en las que el indicador \(MSFT\) tom? el valor m?s alto y los indicadores \(DD\) y \(MCD\) el valor m?s bajo, los dem?s valores son bastante similares a los del segundo cl?ster.
Aplicamos el m?todo de K-medias para las series de tiempo del Dow Jones con tres cl?sters,
grupos<-kmeans(DJTable[,-1],3,iter.max = 1000)
# grupos$cluster
# grupos$centers
par(mfrow=c(3,1))
plot(grupos$centers[1,],type="o")
plot(grupos$centers[2,],type="o")
plot(grupos$centers[3,],type="o")
Al igual que en el ejercicio anterior, las gr?ficas parecen ser bastante similares, por lo que se mostrar?n en el mismo plot de modo de poder compararlas de forma m?s objetiva
par(mfrow=c(1,1))
plot(centros[1,],type="o")
lines(centros[2,],type="o", col=2)
lines(centros[3,],type="o", col=3)
El an?lisis en este caso es exactamente igual al hecho en el ejercicio anterior ya que los grupos identificados al utlizar k-means fueron los mismos que utilizando clasificaci?n jer?rquica. Esto es, el primer cl?ster est? formado por aquellas fechas en las que los indicadores \(CAT\), \(CVX\), \(IBM\), \(KO\), \(MCD\) y \(MMM\) tomaron valores considerablemente m?s altos. El segundo cl?ster est? formado por aquellas fechas en las que los indicadores \(AA\), \(PFE\), \(UNH\), \(VZ\), \(WMT\) y \(XOM\) alcanzaron los valores m?s bajos, por ?ltimo, el tercer cl?ster est? formado por aquellas fechas en las que el indicador \(MSFT\) tom? el valor m?s alto y los indicadores \(DD\) y \(MCD\) el valor m?s bajo, los dem?s valores son bastante similares a los del segundo cl?ster.